Annihilating random walks in one-dimensional disordered media 
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We study diffusion-limited pair annihilation A + A — > on one-dimensional lattices with inhomo- 
geneous nearest neighbour hopping in the limit of infinite reaction rate. We obtain a simple exact 
expression for the particle concentration p^ (t) of the many-particle system in terms of the conditional 
probabilities P(m; t\l;0) for a single random walker in a dual medium. For some disordered sys- 
tems with an initially randomly filled lattice this leads asymptotically to p(t) = P(0; 2t\0; 0) for the 
disorder-averaged particle density. We also obtain interesting exact relations for single-particle con- 
ditional probabilities in random media related by duality, such as random-barrier and random-trap 
systems. For some specific random barrier systems the Smoluchovsky approach to diffusion-limited 
annihilation turns out to fail. 
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Stochastic reaction-diffusion processes play an impor- 
tant role in the description of interacting many-particle 
systems in both physics and chemistry. Usually real sys- 
tems are much too complex to be amenable to analytical 
or even numerical investigation. However, particularly 
in the context of critical phenomena, simple toy models 
may suffice to determine universal properties correctly 
and to predict and explain observed power laws or uni- 
versal amplitude ratios. Hence it is of importance both to 
examine the behaviour of such models and to understand 
possible relationships between microscopically different 
processes and their characterization in terms of univer- 
sality classes. In this letter we investigate a model of 
annihilating random walkers describing diffusion-limited 
pair annihilation (DLPA) of identical particles in inho- 
mogeneous media. This process describes both chemical 
reactions where the particles change their state into an 
inert reaction product which takes no part in the sub- 
sequent dynamics of the system, or physical reactions 
where the particles actually annihilate under the emis- 
sion of radiation. For a recent review of experimental 
and theoretical applications of the ordered model, see 
Rcfs. 0,^]. The model is related by a similarity transfor- 
mation to the diffusion-limited coagulation process ||]4|] 
which describes e.g. laser-induced exciton dynamics on 
polymers. The model also maps to Glauber dynamics for 
the one-dimensional Ising model ^)-@| and is therefore of 
interest in the study of the kinetics of disordered equi- 
librium systems. It plays an important role not only for 
the study of spin-relaxation phenomena but also for the 
solution of more intricate problems such as the derivation 
of persistence exponents ]i| and non-equilibrium steady 
states §,0. 

The first comprehensive treatment of homogeneous 
diffusion-limited annihilation dates back to Smolu- 
chovsky's classical work in 1917 p| , but despite renewed 
strong interest in the 80's and 90's the experimentally 
more realistic case of spatially inhomogeneous particle 
hopping rates, e.g. in the presence of quenched disor- 



der, has so far received very little attention (l|jl|,(§. In 
homogeneous, translationally invariant environments the 
particle density decays with a power law which depends 
on the dimensionality of the system. Both theoretically 
and experimentally one finds p(t) ~ 1 /v ' Dt in one dimen- 
sion, in agreement with exact results Interest- 
ingly, this result is at variance with the (dimensionality- 
independent) mean-field behaviour p(t) ~ l/{Dt) which 
is correct only in three (and higher) dimensions. The am- 
plitude is universal in the sense that it depends neither 
on the initial density (for random initial conditions) nor 
on the reaction rate |lq| . However, if particles are mov- 
ing in an non-translationally invariant energy landscape 
it is not obvious how this will change the decay of the 
local or overall particle concentration. The physical mo- 
tivation behind the study of the one-dimensional case is 
not only its experimental relevance for polymer physics, 
but also its theoretical importance in the understanding 
of the role of fluctuations in low-dimensional systems. In 
both one- and two-dimensional systems diffusive mixing 
is inefficient and leads to the building up of large-scale 
correlations. Thus the classical mean-field rate equations 
for the study of these systems tend to fail and require a 
more sophisticated treatment. 

Before going into details, we remind the reader of the 
essentials of Smoluchovsky's reasoning. The idea behind 
this approach JllJ is to replace in the mean-field rate 
equation for the density p(t) = —\p 2 (t) the reaction 
constant A by an effective time-dependent reaction rate 
which is proportional to the diffusive current j(t) into 
an absorbing particle in a background of constant den- 
sity. In one dimension where j(t) oc yJD/t this leads (up 
to the universal amplitude which cannot be determined 
from the Smoluchovski argument) to p(t) oc l/\/Dt. One 
would like to know whether the fluctuation-improved 
mean-field theory of Smoluchovsky which predicts the 
correct behaviour of the one-dimensional ordered system 
remains valid in the presence of disorder. Naively, one 
might expect that the diffusion constant of the pure sys- 
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tern would have to be replaced by some effective diffusion 
constant of the disordered system, i.e. p(t) oc l/-y/D e fft 
for the disorder averaged density. But clearly this can- 
not always work as can be seen in a simple and natural 
example: Suppose one investigates DLPA on an ensem- 
ble of ordered chains of varying length L < A, modeling 
e.g. a polymer mixture of polymer fragments of vary- 
ing finite length. This is equivalent to taking an infinite 
chain and place randomly, but with maximal distance 
A, broken bonds across which particles cannot move. If 
initially particles are placed randomly with probability 
1/2 on each lattice site, then the steady state density in 
each chain of length L is equal to p* = 1/(2L). Assum- 
ing that each chain length occurs with equal probability, 
then the averaged density p* = 1/(2L) ~ lnA/(2A) for 
large maximal length A. On the other hand, since the 
particle is confined to a box of length L the system is 
subdiffusive with -D e // oc A 2 jt. Hence the Smoluchovsky 
formula gives the wrong result p* oc 1/A. Even worse, 
with this Smoluchovsky approach one cannot even esti- 
mate the approach of the density to its stationary value. 
One could, of course, try to be smart and apply Smolu- 
chovsky's approach to a finite system and then average 
over system size. This gives indeed p* oc lnA/(A). How- 
ever, if applied to a different system where one distributes 
infinitely deep traps (sites out of which particles cannot 
jump) at a maximal distance A, then this refined Smolu- 
chovsky argument would give p* oc In A/ (A) also for this 
model. However, the exact result given below shows that 
in fact for this type of disorder p* = 1/A which hap- 
pens to be consistent with the first, naive Smoluchovsky 
result. We conclude that there is no simple argument 
which tells one how to proceed without already knowing 
the answer. Thus, exact results are required for a study 
of DLPA in inhomogeneous media. 

To this end we investigate by a new exact mapping a 
DLPA lattice model with space-dependent hopping rates. 
The particles have no attractive or repulsive interaction 
between themselves, they hop with fixed rates r k (if.) 
from lattice site k to site k+1 (k — 1). When two particles 
meet on site k they both annihilate instantaneously. This 
limit of infinite reaction rate corresponds to the renormal- 
ization group fixed point of the ordered system J2|,[l6],[l7j 
and we believe that also the disordered system with finite 
reaction rate will be in the same universality class as the 
infinite rate limit. 

In the case of constant hopping rates the predictions of 
this model are in excellent agreement with experimental 
data on exciton dynamics on very long ordered polymer 
chains. Thus we expect that the disordered model gives 
an equally well description of the behaviour of realistic, 
disordered systems. 

We define the process in terms of a master equation for 
the probability P(n] t) of finding, at time t, a configura- 
tion r] of particles on a lattice of L sites. Using standard 
techniques J18|— p0[ we express the time evolution given by 



the master equation in terms of a quantum Hamiltonian 
H. Since particles annihilate instantaneously when they 
meet, there can never be more than one particle on any 
given site. The idea is now to extend the bosonic Fock 
space formalism of ]l8| p0| to a Quantum spin chains rep- 
resentation [|T] ^3). One represents each of the 2 L pos- 
sible particle configurations 77 by a vector 1 77 ) which to- 
gether with the transposed vectors ( r\ | form an orthonor- 
mal basis of a vector space X = (C 2 )® L . A state r] with 
iV" particles placed on sites ki, . . . , fcjv is represented by 
I k±, ■ ■ ■ , A; at ), the completely empty lattice by the vec- 
tor I ) . The probability distribution corresponds to a 
state vector | P(t) ) = J2r/ex P( 7 T>t)\ 7 l) an d one writes 
the master equation in the form 

±P(r ) ;t) = -{r 1 \H\P(t)) (1) 

where the off-diagonal matrix elements of H are the (neg- 
ative) transition rates between states and the diagonal 
entries are the inverse of the exponentially distributed 
life times of the states. A distribution at time t is given 
in terms of an initial state at time t = by | P(t) ) = 
e -ff*|p(0)). The expectation value p k (t) = (s\n k \P(t)) 
for the density at site k is given by the projection oper- 
ator n k which has value 1 if there is a particle at site k 
and otherwise. The vector ( s \ — X^ex ( ^ I performs 
the average over all possible final states of the stochastic 
time evolution. Choosing the basis of X such that a par- 
ticle (vacancy) on site k corresponds to spin up (down) 
the Hamiltonian 

H = -J2(^h++ikK) (2) 
k 

for the process can be written in terms of Pauli matrices 
h k = ( s k s k±i + s k s t±i- n k) wner e n k = (l-crf)/2 and 
s k = (a k ± io v k )j2 create (s k ) and annihilate (s k ) par- 
ticles (see jl| for details). Since the time evolution con- 
serves particle number modulo 2, it is convenient to work 
only on the even and odd subspaces defined by the pro- 
jector P ± = (1 ± Q)/2 where Q = (-l) w = ]J k a z k . For 
averaging over final states we then use (s* | = (s IP^. 
The projection on the even sector of the uncorrelated ini- 
tial state with a density 1/2 used below is given by the 
vector I 1/2+) = (l/2) L - 1 \s+). 

In one dimension with homogeneous nearest neighbour 
hopping DLPA is related to zero-temperature Glauber 
dynamics by a domain-wall duality transformation || 
which is an invertible similarity transformation 0. On 
the other hand, zero-temperature Glauber dynamics can 
be brought by another similarity transformation into a 
form which is the transpose of the Hamiltonian for (ho- 
mogeneous) DLPA (^4|. We use these results to con- 
struct a new matrix V such that the process defined by 
H = V- 1 H T V also describes a DLPA process with near- 
est neighbour hopping albeit with different hopping rates 
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in the presence of disorder. We find for the even particle 
sector 

T>+ = 7172 ■ • ■ 721,-1 (3) 

where j 2 k-i = [(1 + i)cr% - (1 - i)] /2, j 2 k = 
[(1 + i)(T x k a x k+1 - (1 - i)] /2 and X>_ = -£>+cr£ for the 
odd particle sector. To see what happens under the map- 
ping we note that T>± is unitary and transforms Pauli 
matrices as follows: 

^M + ^={ti til <*> 

In the even sector one now finds 

H = -J2(rkK+£ k hi_ 1 ). (6) 

k 

This process is of the same form as the original process 
(|J), but with dual hopping rates Ik = fk , fk— &k+i- 
We shall refer to the environment defined by the dual 
rates as to the dual environment ]2j| . 

In order to make use of the mapping one needs to 
know how a given initial distribution and the observables 
change under the transformation. For the transforma- 
tion laws for states one needs (^), (0) and T>+ 1 \s + ) = 
-i(i-l) L - 1 \0),V + \s+) =i(-i-lY'- 1 \0) for the even 
sector and analogous relations for the odd sector. For 
the density at site k for an arbitrary initial state one 
then finds 

( s+ \n k e- m \ P > = ~< P \Ve-"\\ a%_ x cj%)V- l \ s+). 

(7) 

The transformed initial state is a superposition of the 
steady state (the empty lattice) and the two-particle 
state with particles at sites k — l,k. Bearing in mind 
that H does not have any particle creation terms one 
now realizes that the time-dependence of the density for 
an arbitrary many-particle initial distribution is com- 
pletely given by the dynamics of just two annihilat- 
ing random walkers in the dual disordered environment, 
Pk(t) = a - J2 m ,i a mi(m, I \e~ m \ k-l,k) where the co- 
efficients a, ai m are determined by the initial state and 
straightforward to work out J2(| • 

In order to avoid immaterial technical complications 
with boundary terms, we consider from now on only in- 
finite systems. By choosing some of the hopping rates 
equal to zero one can always recover results for finite 
systems. To calculate the two-particle transition proba- 
bility P(m, n; t\k, I; 0) = (m, n\e~ Ht \ k,l) we note that 
the transition probability for a single random walker 
P(m,k;t) = P(m; t\k,0) is the sum over all paths lead- 
ing from k to m, each weighted with its proper statistical 



weight given by the hopping rates and the particular form 
of the trajectory. Hence, for two non-interacting parti- 
cles moving from k to m and from / to n respectively, 
P(m, n; t\k, I; 0) = P(m, k; t)P(n, I; t). This sum includes 
the contribution of paths which cross each other. In an 
annihilating random walk of otherwise non-interacting 
particles the contribution of all crossing paths have to 
be subtracted. Since we are on a one-dimensional lattice 
this contribution is just the one given by all paths which 
start at site k and end at site n (instead of m) and which 
start at site I and end at site m (instead of n) . Therefore 

P(m, n; t\k, I; 0) = P(m, k; t)P(n, I; t) - P(n, k; t)P(m, I; t) 

(8) 

This further reduces the calculation of the density to the 
solution of a single-particle random walk problem in the 
dual random environment. For an uncorrelated random 
initial state with density 1/2 in the sector of even particle 
number one gets 

p fc (t) = < 2 |e- A '| fc - 1, fe >/2 (9) 

where ( 2 | = ^ n>m (m,n\ is the sum over all states 
with two particles. Thus the density at site k is equal 
to one half the survival probability of two annihilating 
random walkers starting at sites k— 1, k and moving in the 
dual environment. As a first specific result we calculate 
the final density of a system with infinitely deep traps 
placed randomly as discussed above. Clearly, if k is not 
a trap site, then pfc(oo) = 0. Hence the disorder-averaged 
density p(oo) = qp(co) where q — 2/A is the density of 
traps and p(t) is the density at a trap site. From (||) one 
finds p(t) — 1/2 and therefore p(oo) = 1/A. 

Consider now the relation between the single-particle 
conditional probabilities for dual environments. One may 
write P(m;t\k;0) = (s {""ifi^ 1 afn m e- m \ k). By 
taking the transpose and transforming under 2?_ one gets 
in the infinite volume limit the interesting exact relation 

P(m,k;t) - P(m,k- l;t) = 
P(m- l,k- l;t) - P(m - l,k;t). (10) 

Note that eqs. (0) - © hold for any fixed hopping envi- 
ronment, disordered or inhomogencous, but regular. In 
this letter we focus on disordered systems with transla- 
tionally and reflection invariant hopping rate distribu- 
tions. In this case ( |l0| ) gives 

P(^t)=P(r;t) + c(t). (11) 

with an undetermined function c(t) which is irrelevant for 
what follows. Relations (JTc|) , ( |TT| ) are remarkable in that 
they relate the conditional probabilities for dual systems. 

In order to analyze (^) further we take a mean 
field approach to the disorder average, i.e. we re- 
place in (pf) the disorder average of the conditional 
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probability for two distinguishable non-interacting ran- 
dom walkers P(m, k; t)P(n, /; t) by the factorized average 

P{m,k;t) P(n,l;t). 

We have convinced numerically that this factorization 
holds well for the random barrier model with uncor- 
rected bond hopping probabilities bk drawn uniformly 
from the interval 0.05 < bk < 1/2. This was done by ex- 
act numerical solution of the discrete-time master equa- 
tion for a random walker on a lattice of L = 200 sites 
for a given random realization of the disorder and then 
taking the average over 100000 disorder realizations (to 
keep disorder-related fluctuations small). To show this 
factorization in Fig. 1 the function 



Ki k (t) 



P(m,k;t)P(n,l;t) 
P(m,k;t) P(n,l;t) 



(12) 



is shown for m = n = 1 and k = I = in a double 
logarithmic plot. The return probability gives the largest 
contribution to the function R™ k (t), hence we do not 
need to check the factorization at other positions. The 
function R^it) fitted well by a power law 
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where a = —0.49. For this the fit routine of the soft- 
ware package Mathematica was used. We expect that 
the exact exponent to be a = 1/2. The deviation of the 
last point in Fig. 1 from the line originates from finite- 
size effects. The square root of the time at this point 
is about the length of the system and therefor the con- 
stant stationary probability of the random barrier system 
is reached exponentially fast. We have also investigated 
the case where zero is the lower limit of the distribution 
of the hopping probabilities in the same system as above. 
In this case the behaviour of the mean square displace- 
ment of a single particle is sub-diffusive < x 2 >~ tj log(t) 
p7| , |30"[ . This is due to the existence of very small hopping 
probabilities in even very small regions in the chain. The 
factorization can only hold for very long times. Never- 
theless the function R™ k {t) is very small even for short 
times and for longer times there is a tendency of the func- 
tion R™ k {t) to decrease as it can be seen in Fig 2 where 
R™ k (t) is spread vs. log(t). Independent arguments for 
the validity of this assumption for more general types of 
disorder are given in the conclusions. 

Then with (||) and ( |i"T| ) the disorder average of the 
density for the initially randomly filled lattice is given by 



p(f) = (p(0;2i)+P(l;2i)) /2 



(14) 



For large times (|14|) becomes the return probability 
quoted in the abstract. Having in mind processes like 
exciton dynamics it is reasonable to consider (A) ran- 
dom bond disorder ru = £k+i = bk where hopping across 
a bond k, k + 1 is symmetric, but bond-dependent and 



(B) random site disorder rk = £k = Sfc- The energy of 
a particle in the random barrier model (A) is the same 
at each site, but between sites there are energy barriers 
of random height Ek- Thermal fluctuations cause the 
particle to jump over these barriers with a random rate 
bk oc exp (— f3Ek). Case B corresponds to the random 
trap model. Here the particle sits in a site-dependent 
potential of depth — Ek- Random bond and random site 
disorder are dual in the sense of Eq. (^). Since for 
the random-bond model P(m, k; t) = P{k, m; i) we con- 
clude that the disorder-averaged conditional probabilities 
of the random-barrier model and the random-trap model 
are equal in one dimension for any translationally in- 
variant disorder distribution up to a function c(t). For 
an uncorrelated ergodic disorder distribution, diffusion in 
random barrier systems converges to Brownian motion, 
i.e. the averaged conditional probability becomes asymp- 
totically equal to a Gauss distribution with an effective 
diffusion constant 



(15) 



D eff=K 

and therefore asymptotically 



P(t) 
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(16) 



for random trap and random barrier systems. For ex- 
ponentially distributed barrier or trap energies v(E) oc 
exp (— E I a) the random walk becomes subdiffusive be- 
low a critical temperature given by a(3 c = 1 pl| , ^2[ . 
This leads to a time-dependent effective diffusion co- 
efficient D e fj oc t( 1 ~ cr ' 3 )/( 1 + <T ' 3 ) and hence to a slower 
non-universal power law decay of the density pit) oc 

t -l/(l+<T(9)_ 

To conclude, we obtained the following new results. 

(i) We found the duality relations ([loj), ([ll]) for ran- 
dom walkers in dual systems. They play a role in the 
derivation of (jlj) - (|lq), but are also of interest in their 
own right. 

(ii) The Smoluchovsky approach to DLPA is consistent 
with random trap and random bond systems with van- 
ishing broken bond probability ([l5]), (|l6|), but fails for 
some important disordered broken bond systems in one 
dimension. This raises the question under which general 
conditions on the disorder the Smoluchovsky approach is 
correct (see below). 

(iii) The expectation value of the density in diffusion- 
limited annihilation in one dimension with nearest neigh- 
bour hopping and infinite annihilation rate is completely 
determined by the dynamics of a single random walker 

(H) m a dual hopping environment (^J). This exact 
result allows for the exact calculation of the density in 
specific environments, but can also be used for extremely 
accurate numerical calculations of the density for any 
fixed inhomogeneous hopping environment or for disor- 
dered systems with subdiffusive logarithmic behaviour. 
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The approach has a straightforward extension to the cal- 
culation of correlation functions. 

(iv) For disordered environments which are on average 
translationally and reflection invariant and which lead 
to an asymptotic factorization of the two-particle condi- 
tional probabilities, the density at time t with a random 
initial state is equal to the return probability ( |l4| ) of a 
single particle. An important open problem is the deriva- 
tion of conditions on disorder distributions under which 
factorization holds. For such distributions (|l4|) is an ex- 
act asymptotic result, and, as we would like to point out, 
consistent with the Smoluchovsky approach if the return 
probability is proportional to 1/ U D e fft. This may be a 
hint under which circumstances the Smoluchovsky treat- 
ment is adequate for the calculation of the density. Since 
at late times particles are separated (on average) by a 
distance 1/p — > oo, one would not expect a finite reac- 
tion rate or short-range interactions between particles to 
change this asymptotic behaviour. It is interesting to 
note that (^) is also consistent with the renormalization 
group treatment of DLPA with disorder (u|] . Therefore 
we are confident that (Q) holds for a broad class of dis- 
order distributions. 
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Figure Captions 



FIG. 1. Factorization of one particle condition- 
probabilities after the disorder average. R(t) vs. t is 
shown in a double logarithmic plot. The line has the slope 
—0.49 and is a fit of the first 7 points which are produced 
by exact numerical calculations of the master equation 
where the jump-probabilities bk were taken randomly 
from a uniform distribution between 0.05 < bk < 1/2. 

FIG. 2. Factorization of one particle condition- 
probabilities after the disorder average. R(t) vs. log(t). 
The points are produced as in Fig 1, but the jump- 
probabilities bk were taken randomly from a uniform dis- 
tribution between < bk < 1/2. 
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